Topological Speed Limits to Network Synchronization 



m 
o 
o 

(N 

C 
3 



Marc Timme, Fred Wolf, and Theo Geisel 
Max- Planck- Institut fur Strdmungsforschung, 37073 Gottingen, Germany 

We study collective synchronization of pulse-coupled oscillators interacting on asymmetric random 
networks. We demonstrate that random matrix theory can be used to accurately predict the speed 
of synchronization in such networks in dependence on the dynamical and network parameters. 
Furthermore, we show that the speed of synchronization is limited by the network connectivity and 
stays finite, even if the coupling strength becomes infinite. In addition, our results indicate that 
synchrony is robust under structural perturbations of the network dynamics. 
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Complex networks have attracted an enormous re- 
search interest in the recent past pj. Most studies have 
focussed on the structures of a variety of systems such 
as the world wide web, email networks, genetic networks, 
and biological neural networks 0, An equally im- 
portant task is to understand the collective dynamics on 
such networks. In particular, the question arises: how 
is the dynamics on a complex network influenced by its 
structure Q? 

Synchronization appears to be one of the simplest 
kinds of collective dynamics among coupled dynamical 
systems 0. Q|. It was found to be ubiquitous in arti- 
ficial as well as natural networks as different as Joseph- 
son junction arrays and biological neural networks]^. 
To understand the dynamics of such networks, theoreti- 
cal studies have emphasized systems consisting of simple 
units such as phase- and pulse-coupled limit-cycle oscilla- 
tors Yet, although most real- world networks often 
possess a complex connectivity structure, most studies 
of synchronization of coupled oscillators are either re- 
stricted to networks of globally coupled units and sim- 
ple regular networks, or work in some mean field limit 
[3, l| • Although exact results on synchronization in net- 
works with a general structure have been obtained re- 
cently (slbiillll; it is still not well understood how the 
structure of a complex network affects dynamical features 
of synchronization. 

In this Letter, we study the collective synchronization 
of pulse-coupled oscillators interacting on asymmetric 
random networks. We find that the speed of synchroniza- 
tion is restricted by the network connectivity and stays 
finite, even if the coupling strength becomes infinite. No 
such speed limit exists in large networks of globally cou- 
pled units. More generally, we show that the theory of 
random matrices can be used to successfully predict the 
speed of synchronization as a function of dynamical and 
network parameters. In addition, our results indicate 
that synchrony occurs robustly, i.e. persists under struc- 
tural perturbation of the network dynamics. 

We consider asymmetric random networks of N oscil- 
lators which interact by sending and receiving pulses [l^ . 
The sets Pre(z) of presynaptic oscillators that send pulses 



to oscillator i specify the structure of such a network. For 
each oscillator i, the fc; := | Pre(i)| presynaptic oscillators 
are drawn from the uniform distribution among all other 
oscillators {1, . . . , A r }\{i}. 

A phase- like variable 4>i(t) specifies the state of each 
oscillator i at time t. In the absence of interactions, 
the dynamics of unit i is given by dcpi/dt = 1. When 
oscillator i reaches a threshold, <j>i(t) = 1, its phase 
is reset to zero, <j>i(t + ) = 0, and the oscillator is said 
to 'fire'. A pulse is sent to all postsynaptic oscilla- 
tors j G Post(i) which receive this signal after a de- 
lay time r. The incoming signal induces a phase jump 
<f>j({t+r) + ) := U~ 1 {U((f) : j(t + T))+£ j i) which depends on 
the instantaneous phase (t + r) of the postsynaptic os- 
cillator and the coupling strength Eji which we take to be 
inhibitory (phase- retarding), Sji < 0. The phase depen- 
dence is determined by a twice continuously differentiable 
'potential' function U((j>) that is assumed to be strictly 
increasing, U'(<p) > 0, concave (down), U"{4>) < 0, and 
normalized such that U(0) = 0, U(l) = 1 (cf. 0). We 
focus on the specific form U(4>) = Uif(4>) = 1(1 — e~ TlF ^) 
that represents the integrate-and-fire oscillator defined 
by the differential equation V = I — V (and a thresh- 
old at V = 1). Here I > 1 is an external input and 
Tif = log(// (I — 1)) the intrinsic period of an oscillator. 
Other forms of U(<fi) give qualitatively similar results. 
In such a network the synchronous state, <f>i(t) = <f>o(t) 
for all i, exists if the coupling strengths are normal- 
ized such that YljeFie(i) £i i = £ - ^ s P er i°d is given by 
T = t + 1- U- 1 (U(t)+s). 

In numerical simulations of the network dynamics, we 
find that the synchronous state is always stable, indepen- 
dent of the parameters (cf. P3|)- A sufficiently small per- 
turbation 6(0) = d = (Si, . . . , Sn) t of the phases, defined 
by <j>i (0) = <^>o(0)+£i asymptotically decays exponentially 
with time. Thus, denoting 8'(t) := S(t) — lim s _voo <5(s), 
the distance A(n) := max^ |<5-(nT)|/ max^ |<5-(0)| from the 
synchronous state behaves as 



A(n) ~ exp(-n/r syn ) 



(1) 



defining a synchronization time r syn in units of the collec- 
tive period T. The speed of synchronization r™ strongly 
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depends on the parameters. For instance, as might be 
expected, synchronization is faster for stronger coupling. 
Surprisingly, however, we find that synchronization can- 
not be faster than an upper bound even if the coupling 
strength becomes arbitrarily large (cf. Fig. . 
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FIG. 1: Asymptotic synchronization time in a random net- 
work (N = 1024, ki = k = 32, I = 1.1, t = 0.05, e tj = e/k 
for j e Post(i)) . The inset shows the distance A of a per- 
turbation S from the synchronous state versus the number of 
periods n (s = —0.4). Its slope yields the synchronization 
time T syn shown in the main panel as a function of coupling 
strength |e|. Simulation data (O)i theoretical prediction ( — 
— ) derived in this Letter, its infinite coupling strength asymp- 
tote ( ). 



representing a uniform phase-shift and thus reflecting 
time-translation invariance. Furthermore, the Gersh- 
gorin theorem 0] implies that all eigenvalues are located 
inside a disk of radius tq = 1 — A centered at A , such 
that in particular |A»| < 1 and the synchronous state is 
(at least neutrally) stable. For simplicity, we consider 
networks of homogenous random connectivity, ki = k for 
all i, in the following. 

We numerically determined the eigenvalues of differ- 
ent stability matrices A for various network sizes N € 
{2 6 , . . . , 2 14 }, in-degrees k e {2, . . . 2 8 }, and dynamical 
parameters e, t, and I such that < Ao < 1. In general, 
we find that for sufficiently large k and N the non-trivial 
eigenvalues resemble a disk in the complex plane that is 
centered at about Aq and has a radius r that is smaller 
than the upper bound given by the Gershgorin theorem, 
r < 1 — Aq. Examples are shown in Fig. [21 
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To understand how the speed of synchronization de- 
pends on the dynamical and network parameters, we 
analyze the linear stability of the synchronous state. 
Following [jjj we obtain a nonlinear stroboscopic map 
S(T) — F(S) for the perturbations, the linearization of 
which reads 



S(T) = AS 



(2) 



where the elements of the stability matrix A are given 
by Aij = p it „ - p i<n -i if j = j n G Pre(i), An = 
Pifi, and A^ = otherwise. Here j n identifies the 
n th pulse received during this cycle by oscillator i and 
p i>n := U'(U-\U(t) + EZ lr=l e ijm ))/U'(U- 1 (U(T)+e)) 
for n G {0, 1, ... , h}. For U{<j>) = Uiy{4>) and coupling 
strengths e^a = e/ki for j £ Pre(«) the matrix elements 
reduce to [13] 



Ai 



l ~ Ao if j G Pre(i) 
if j = i 

if j i Pre(i) U {i} 




where 



A 



Ie 



-tT u 



Ie 



~tT, 



IF _ £ 



> 0. 



(3) 



(4) 



Obviously, A constitutes a row-stochastic matrix, i.e. 
J2j Aij — 1 for all i. Thus A has one trivial eigenvalue 
Ai = 1 associated with the eigenvector V\ = (1,1,..., 1) T 



Figure 2: Distribution of eigenvalues Ai of two stability ma- 
trices A in the complex plane (/ = 1.1, e = —0.2, r = 0.05 
^A w 0.83; k = 8) for networks of (a) N = 32, (b) N = 512 
oscillators. For large networks, the non-trivial eigenvalues 
seem to be distributed uniformly on a disk in the complex 
plane. The prediction from random matrix theory (Eq. {S3} ) 
is indicated by a circle. The arc through the trivial eigenvalue 
Ai = 1 is a sector of the unit circle. 



This eigenvalue distribution is reminiscent of the "cir- 
cle law" of random matrix theory 0|: Gaussian asym- 
metric random matrices, having a distribution of matrix 
elements 



PGauss 



2 exp 



N4 
2r 2 



(5) 



with independent Jij and Jji , also exhibit eigenvalue dis- 
tributions 



PGauss 

(A) 



(Trr 2 )- 1 




if |A| < r 
otherwise 



(6) 



for N — > oo that are uniform in a disk in the complex 
plane ^| . The radius r of the disk is given by 



(7) 



where a 2 — (JB) is the variance of the matrix elements. 

Interestingly, we find that the radii of the eigenvalue 
distributions of the above stability matrices (j^J well 
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agree with the radii obtained from Eq. if (Jfj) is 
replaced by the variance of the elements of the stabil- 
ity matrices shifted such that they also exhibit a zero 
average eigenvalue. To directly compare the eigenvalues 
of the stability matrices, which have average eigenvalue 
[Aj] := jr Yli=i = A o to those of the Gaussian ensem- 
ble, we transform A!^ — Aij — S^Aq shifting the average 
eigenvalues to [A-] = 0. Here Sij denotes the Kronecker 
delta, Sij = 1 if i = j and Sij — if i ^ j. For the 
variance of A' we obtain 



[ A 'tj 2 ] ~ [ A 'ij] 



1 

N 



N 

£4 



N 



(8) 



(9) 



For identical non-zero coupling strengths, the off- 
diagonal sum is exactly equal to J2j=i j^i A 2 j = (1 — 
Ao) 2 /k such that, using 0, we obtain the random ma- 
trix theory prediction 



fRMT = Ni(J A , = (1-A )[-- — 



(10) 



for the radius r of the disk of eigenvalues of the stability 
matrices A 

We verified this scaling law ifTTijl for various dynami- 
cal parameters Aq (determined by different /, e, and r), 
network sizes N, and in-degrees k and found excellent 
agreement with numerically determined eigenvalue dis- 
tributions, see, e.g., Fig. To quantify the accuracy of 
the prediction (fflTlfl we numerically estimated the radius 
of the distribution of the non-trivial eigenvalues of A for 
various N, k as well as Ao. Results from two different 
estimators are shown in Fig. 03 The real part estima- 
tor rR e := | (max^i Re(Ai) — min^i Re(Ai)) estimates 
the radius from the maximum spread of eigenvalues par- 
allel to the real axis. Typically, r^ e should give an es- 
timate that is slightly inaccurate because it is based on 
two eigenvalues only. This is circumvented by the average 
estimator r av := |A; - {A - (1 - A Q )N- 1 )\ 

that estimates the radius r of a circle from the average 
distance (d) of eigenvalues from its center, because (d) = 
Jo* fo r ' 2 p( r ')drd<p = |r if we assume a uniform density 
p(r') according to Q. Here we take the center of the disk 
to be the average (A;} i>2 = A - (1 - A )N- 1 + 0(N~ 2 ) 
of the non- trivial eigenvalues. Varying k at fixed TV as 
well as N at fixed k yields excellent agreement between 
the numerical data and the theoretical predictions for 
sufficiently large N and k (Fig. . Varying the coupling 
strength |e| and thus Aq yields equally good agreement 
(cf. Fig.[TJ). 

The radius iflOj) implies a prediction for the synchro- 
nization time (see (QJ) 
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Figure 3: Scaling of the radius r of the disk of non-trivial 
eigenvalues. Main panel displays the radius r as a function 
of network size N for fixed k = 32. Symbols display and rR e 
(x) and r av (CO- Inset displays r as a function of k for fixed 
N = 1024. Dots display numerical data of r av . In the main 
panel and the inset, lines are the theoretical estimate trmt 
(Eq. IWt ) Other parameters as in Fig. [5] 



in terms of the (in modulus) largest non-trivial eigen- 
value A m . With increasing coupling strength the 
synchronization time decreases. However, the speed of 
synchronization r a y* is bounded by a finite speed for ar- 
bitrary large |e|: Even if |e| ^> 1 and thus Aq < 1, 
the largest non-trivial eigenvalue asymptotically becomes 
A m « fc -1 / 2 for sufficiently large N. Thus the shortest 
synchronization time 



syn 



2 

LTfc 



(12) 



'syn 



-l/ln(A + r RMT ) , 



(11) 



is limited by the network connectivity (cf. the asymp- 
tote in Fig. nj. This means that even for arbitrary 
strong interactions, the speed of synchronization stays 
finite. Furthermore, at fixed k, the synchronization time 
also cannot exceed a certain maximum, even if the net- 
work size N becomes extremely large (cf . Fig. . This 
bound t^^*°° is determined by the asymptotic radius 
^ := liniAT^oo r RMT = (1 - ylo)fc~ 1/2 . Moreover, be- 
cause eigenvalues change continuously with a structural 
perturbation to the system's dynamics, the existence of 
a gap g := 1 — (Ao + ^00) > indicates that no eigenvalue 
crosses the unit circle for sufficiently small structural per- 
turbations. Thus stable synchrony is not restricted to the 
specific model considered here, but persists in systems 
obtained by structural perturbations of the dynamics. 

The above results show, that the distribution of eigen- 
values of a sparse stability matrix with deterministic 
non-zero entries at certain random positions is well de- 
scribed by the eigenvalue distribution of the Gaussian 
ensemble, which consists of fully occupied matrices with 
purely random entries. This sparse-Gaussian coincidence 
for asymmetric matrices is similar to that of symmetric 
random matrices for large k: Gaussian symmetric ma- 
trices exhibit an eigenvalue distribution Pa auss W, t ne 
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Wigner semicircle law [2flj|. Sparse symmetric matrices 
exhibit an eigenvalue distribution Ps parse (A) that is 
different from the semicircle law but approaches it in the 
limit k — » oo. For k 3> 1 the distribution of eigenval- 
ues of sparse asymmetric random matrices Pg parse appear 
to be well approximated by the eigenvalue distribution 
of Gaussian asymmetric matrices, Ps parse (^) ~ PGaussM- 
Our results indicate, that this is true even for moderate 
k w 10. 

Further investigations of eigenvalue distributions for 
small- world networks show that with decreasing random- 
ness the speed of synchronization decreases (the second 
largest non-trivial eigenvalue increases) such that com- 
pletely random networks synchronize faster than small- 
world networks, at least asymptotically. 

In conclusion, we have derived accurate analytical pre- 
dictions for the (asymptotic) speed of synchronization 
in asymmetric random networks of oscillators in depen- 
dence of the dynamical parameters s, t, I, as well as 
the network parameters N and k. Even the scaling with 
network size TV, artificially introduced via the variance 
© of finite matrices, is also accurately reproduced (see 
e.g. Fig. I2J. As a particular application, we explained 
the intriguing phenomenon, that the speed of synchro- 
nization stays finite, even if the coupling strengths be- 
come infinite. It turned out that the speed is restricted 
by the network connectivity. In addition, at fixed param- 
eters the synchronization time does not increase above a 
finite threshold if the network becomes very large. Fur- 
thermore, the existence of a gap between the non-trivial 
eigenvalues and the unit circle indicates that stable syn- 
chrony is a robust form of collective dynamics in a large 
class of systems. 

Random matrix theory has previously been applied to 
various physical systems that exhibit certain symmetries 
- such as time-reversal symmetry - but an otherwise un- 
known structure. For instance, correlations of energy lev- 
els in nuclear physics and quantum mechanical properties 
of classically chaotic systems have been successfully pre- 
dicted (see Ref. 0] for a recent review). Our results 
demonstrate that random matrix theory also is an ap- 
propriate tool for analyzing synchronization in random 
networks of dynamical units. Possible lines for future 
applications may include synchronization phenomena of 
pulse- and phase-coupled units as well as of chaotic dy- 
namical systems (cf. also (HI)- The limits of synchro- 
nization speed predicted in this Letter, are expected to 
occur in those systems, too. More generally, other equi- 
libration processes and the dynamics in more structured 
topologies such as small-world networks may be analyt- 
ically investigated using statistical spectral properties of 
the respective operators, too. 

We thank T. Kottos, P. Miiller, H. Sompolinsky, M. 
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